function [dPdv,dRdv] = gradP3(H,impossible,K,Nu,phi,rho,s,tau,P,Vb)

dRdv = gradR(H-1,impossible,K,Nu,phi,rho,tau,P);        % 22 x 10 x Nu for each type
dPdv = cell(1,tau);

for idx1 = 1:tau
                        
    dRdv0 = dRdv{idx1};                         
    
    u = zeros(H,K,Nu);
    dPdv{idx1} = zeros(H+1,H-1,K,Nu);
    
    for idx2 = 1:Nu
        
        dudv = zeros(H,H-1,K);       
        
        for idx3 = 1:K
            
            u(:,idx3,idx2) = (Vb{idx1}(1,idx3,idx2) - Vb{idx1}(2:end,idx3,idx2))./phi(2:end,idx3);
            
            dudv(1:H-1,:,idx3) = repmat(1-impossible(2:H,idx3),[1,H-1]).*(-eye(H-1) ...
                - rho*(repmat(dRdv0(:,1,idx2)',[H-1,1,1])...
                - repmat(dRdv0(:,min(K,idx3+1),idx2)',[H-1,1,1])));
            dudv(H,:,idx3) = -rho*(dRdv0(:,idx3,idx2)' - dRdv0(:,min(K,idx3+1),idx2)');

            for j = 1:H+1

                if impossible(j,idx3) == 1
                        
                    dPdv{idx1}(j,:,idx3,idx2) = grad_dP0dv2(dudv(:,:,idx3),H,...
                        impossible(2:end,idx3),s,u(:,idx3,idx2));
                                            
                else
                        
                    dPdv{idx1}(j,:,idx3,idx2) = integral(@(y)grad_dPdv2(y,...
                        dudv(:,:,idx3),H,impossible(2:end,idx3),j,...
                        phi(2:end,idx3),s,u(:,idx3,idx2)),0,Inf,'ArrayValued',...
                        true,'AbsTol',1e-12,'RelTol',0);                       
                end
            end
            
        end
        
    end
end
